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The Lennard-Jones potential is widely used to describe the interlayer interactions within layered 
materials like graphene. However, it is also widely known that this potential strongly underestimates 
the frictional properties for layered materials. Here we propose to supplement the Lennard-Jones 
potential by a Gaussian-type potential, which enables more accurate calculations of the frictional 
properties of two-dimensional layered materials. Furthermore, the Gaussian potential is computa¬ 
tionally simple as it introduces only one additional potential parameter that is determined by the 
interlayer shear mode in the layered structure. The resulting Lennard-Jones-Gaussian potential is 
applied to compute the interlayer cohesive energy and frictional energy for graphene, M0S2, black 
phosphorus, and their heterostmatures. 
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I. INTRODUCTION 

In 1924, Lennard-Jones published a 12-6 pairwise po¬ 
tential to describe the van der Waals interaction between 
two atoms, which is now known as the Lennard-Jones 
(LJ) potential^ The LJ potential depends only on the 
distance (r) between two interacting atoms. For a rela¬ 
tive displacement u between these two atoms, the varia¬ 
tion in the distance is dr = u • e r , with e r = r/r, which 
shows that this potential is able to effectively control 
the cohesive motion between two atoms. However, the 
LJ potential cannot describe the frictional motion of two 
atoms, because a weak relative frictional motion between 
two atoms does not alter their distance. 

This friction issue can be greatly amplified when the 
LJ potential is applied to the interlayer interaction in 
quasi-two-dimensional layered materials such as bilayer 
graphene^ In these layered structures, the van der 
Waals interlayer interaction is much weaker than the 
covalent intra-layer interaction^ leading to two distinct 
characteristic types of motion in these layered materials, 
i.e., the relative cohesive motion and the frictional mo¬ 
tion. In bilayer graphene, the LJ potential can describe 
the cohesive motion accurately, but it is not able to pro¬ 
vide an accurate measure of the frictional energy^ 

There are only two parameters in the LJ potential - 
one (a) is a length parameter determining the interlayer 
spacing for bilayer graphene, while the other (e) is an 
energy parameter. However, bilayer graphene has two 
independent interlayer motions, i.e., the cohesive motion 
and the frictional motion. As a result, it is not surprising 
that the LJ potential cannot describe both the cohesive 
and frictional motion simultaneously. Several works have 
shown that this friction issue in the LJ potential for bi¬ 
layer graphene can be eliminated by introducing seven 
more potential parameters £ - — 


The aim of the present work is to present a concise sup¬ 
plement for the L J potential in layered materials while in¬ 
troducing a minimum number of fitting parameter, with 
the specific goal of accurately capturing the frictional 
motion. This would be computational beneficial as it 
can be readily implemented in most atomistic simula¬ 
tion packages that use the LJ potential. On the other 
hand, many advanced properties have been found for 
the layered materials, which have garnered both aca¬ 
demic and industrial attention. For instance, few-layer 
graphene can serve as an ideal platform for the investi¬ 
gation of some dimensional crossover phenomena^ - — It 
was found that heterostructures like graphene/MoS 2 can 
mitigate the less desirable properties of each individual 
constituent ] n i 12 Hence, it is important to describe the 
interlayer energy for layered materials more accurately, 
including the important frictional properties^ 

In this paper, we propose to combine the LJ poten¬ 
tial with a Gaussian-type potential (LJ-G) to describe 
the interlayer energy for layered materials. The Gaus¬ 
sian potential introduces only one additional parameter, 
which is determined by the interlayer shear (C) mode in 
the layered structure. The LJ-G potential thus has mini¬ 
mum number of potential parameters and can be applied 
to compute the cohesive energy and frictional energy in 
graphene, MoS 2 , black phosphorus (BP), and their het- 
erostructures. Due to intrinsic lattice mismatch, the fric¬ 
tional energy in all heterostructures is found to be one 
order lower than the individual constituent. 


II. INTERLAYER POTENTIAL 

Fig. [T| shows the AB-stacking for bilayer graphene of 
dimension 30 x 30 A. Both top view and side view are 
shown in the figure. The z-axis is perpendicular to the 
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FIG. 1: (Color online) Top and side views for bilayer graphene 
of dimension 30 x 30 A. 



FIG. 2: (Color online) The energy variation for bilayer 
graphene with respect to different interlayer spacings. dz = 0 
corresponds to the equilibrium interlayer spacing of 3.365 A. 
The interlayer energy is described by the LJ potential with 
parameter values given below Eq. 0- The bottom inset 
shows a zoom in of the curve around the minimum energy 
point, where the energy variation is fit to a quadratic func¬ 
tion E — ^/mz 2 (dz) 2 . /jl = ragra/2 is the effective mass for two 
vibrating graphene layers, with m gr a as the mass of a single¬ 
layer graphene. The fitting parameter u = 88.4 cm -1 gives 
the vibration frequency for the B mode as shown in the top 
inset. 


graphene plane. The x direction is along the armchair 
direction, while the y-direction is in the zigzag direction. 

The following LJ potential is applied to describe the 
interlayer energy, 


V LJ = 4e 


/ <T \ 12 /<T> 

V r / V r > 


(i) 


where r is the distance between two interacting atoms, 
e = 2.96 meV and a = 3.382 A are potential parameters. 
Specifically, the length parameter a is fit to the out-of- 
plane lattice constant in bulk graphite, while the energy 
parameter e is fit to the interlayer breathing (B) mode 
in bilayer graphene, as will be described in further detail 
below. 

To explore the relationship between e and B mode, 
we need to calculate the cohesive energy in the bilayer 
graphene. The structure is optimized via energy min¬ 
imization, after which the cohesive energy for bilayer 
graphene can be computed by evaluating the energy for 
different separation distances of the individual layers. 
Fig. [2] shows the cohesive energy for bilayer graphene. 
The x-axis (dz) is the variation in the interlayer spacing 
with respect to its equilibrium value, so dz = 0 corre¬ 
sponds to the optimized interlayer spacing. 

From lattice dynamical analysis^ it can be shown that 
the strength of the relative cohesive motion is related to 
the frequency of the B phonon mode in bilayer graphene. 
The vibration morphology of the B mode is shown in 


TABLE I: LJG parameters for bilayer graphene. Num¬ 
bers in the parentheses are experimental out-of-plane lattice 
constant and frequency.— The first row is for the LJ po¬ 
tential (g = 0), while the second row is the LJ-G potential. 
Energy parameter is in meV. Length parameter is in A. Fre¬ 
quency is in cm -1 . 


e 

a 

9 

c (6.7) 

c ob (89.5) 

OJC (37.1) 

2.96 

3.382 

0 

6.73 

88.4 

7.3 

2.96 

3.382 

94.87 

6.73 

88.4 

37.1 


the top inset of Fig. [2j More specifically, for the cohe¬ 
sive energy curve around the minimum energy minimum 
(dz = 0), the structure deviates only slightly from its op¬ 
timized configuration. We can thus consider this small 
cohesive motion as a linear vibration of the B mode. 
Hence, we can extract the frequency for the B mode 
from the cohesive energy as shown in the bottom inset 
of Fig. [2j by fitting the cohesive energy to the quadratic 
function E = (1/2 )/iuj 2 (dz) 2 . The quantity fi = ra gra /2 
is the effective mass for two vibrating graphene layers, 
with ra gra as the mass for a single layer of graphene. The 
fitting parameter uj yields the B mode’s frequency. 

We thus fit parameter e in the LJ potential to the fre¬ 
quency of the B mode in the bilayer graphene. Tab. [T] 
shows the fitted LJ parameters for bilayer graphene, 
where the length parameter a in the LJ potential is fit to 
the out-of-plane lattice constant in graphite. The fitted 
LJ potential yields c j = 88.4 cm -1 for the B mode and 
the out-of-plane lattice constant c = 6.73 A. These results 
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FIG. 3: (Color online) The energy variations for bilayer 
graphene as a function of the relative displacement of two 
layers in the x (top) and y (bottom) directions. The inter¬ 
layer interaction is described by the LJ potential. Left insets 
in both panels show the zoom-in of the small displacement 
regime, where the energy variation is fitted to quadratic func¬ 
tions E = 1 fiuj 2 (dx) 2 and E — | (iuj 2 (dy ) 2 . The fitting pa¬ 
rameter uj — 7.3 cm -1 gives the vibrational frequency for the 
C x and C y modes as shown in the right insets. 


are in good agreement with the experimental results i 15 i 16 
It should be noted that the experimental frequency for 
the B mode in bulk graphite (c^buik) has been used to 
extract the frequency of bilayer graphene (c^bi) through 
c^bi = ^buik/V^ according to the linear chain modeler— 
For instance, experiments found c^buik = 126.6 cm -1 in 
graphite^ so the frequency of the B mode in bilayer 
graphene is c^bi = 89.5 cm -1 . This number is listed in 
parentheses in the first line of Tab. U 

Using the above fitted LJ potential, we can also calcu¬ 
late the interlayer frictional energy between two graphene 
layers. Fig. [3] shows the frictional energy for the rel¬ 
ative shearing of two graphene layers along the x and 
y directions. Similar as the cohesive energy, the fric¬ 
tional energy is also in close relation with the interlayer 
phonon modes in bilayer graphene. The frictional en- 
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FIG. 4: (Color online) Interlayer LJ potential for a carbon 
atom on top of a single layer of graphene. The six interlayer 
first-nearest-neighbor van der Walls bonds make most impor¬ 
tant contribution to the interlayer energy. The color bar is 
for the interlayer potential. 


ergy curve around the minimum point determines the 
frequency of the C mode, which is shown in the right top 
inset in both panels of Fig. [3] The frictional energy in 
the x-direction around the minimum point can be fit to 
the quadratic function E = (l/2)/ia; 2 (dx) 2 , in which uo 
is the frequency of the mode (with vibration along 
the x-direction). Similarly, the frictional energy in the 
y-direction gives the frequency of the C y mode (with vi¬ 
bration along the y-direction). For bilayer graphene, we 
find that ujcx = u Cy = 7.3 cm -1 , which is smaller than 
the experimental valued by a factor of |. It implies that 
the frictional energy will be underestimated by a factor 



We have learned that the LJ potential is able to ac¬ 
curately describe the interlayer spacing and the cohesive 
energy between graphene layers. However, this potential 
has a friction issue; i.e., it underestimates the interlayer 
frictional energy in the bilayer graphene by one order. 
This is actually quite reasonable. Considering that there 
are only two parameters in the LJ potential, the predic¬ 
tion of this potential should be limited to two indepen¬ 
dent quantities only, i.e., the cohesive energy (B mode) 
and the interlayer spacing. We should not expect a good 
prediction for the third quantity of interest, the frictional 
energy (C mode). A straightforward solution for this fric¬ 
tion issue is to increase the number of parameters in the 
potential model. For instance, seven additional parame¬ 
ters are introduced in Ref. @ to resolve the friction issue. 

Before presenting our approach, we first make an ex¬ 
plicit examination of this friction issue in the LJ poten¬ 
tial. Fig. [4] shows that for a particular carbon atom 
from the top graphene layer, the LJ potential for this 
atom is mainly contributed by its six interlayer first- 
nearest-neighbor atoms in the bottom layer. We intro¬ 
duce an angle # to describe the direction of these in¬ 
terlayer van der Waals bonds as shown in Fig. SJ For 
bilayer graphene, we have tan# = b/c ~ 0.42, with 
b = 1.42 A as the chemical C-C bond length in the 
graphene plane and c = 3.35 A as the interlayer spac¬ 
ing. We get cos# = 0.92 and sin# = 0.15. For cohesive 
motion, the relative displacement between two graphene 
layers is u = ue z . The resulting variation in the distance 
is dr = ue z • e r = i^cos# = 0.92 u, where e r = f/r is the 
unit vector between two interacting atoms. For frictional 
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FIG. 5: Gaussian shaped potential. Vg = ge p , with p — 
y/x 2 + y 2 as the projection of the distance onto the xy-plane. 
g is the height of the potential. 

motion (eg. in the x-direction), the relative displacement 
between two graphene layers is u = ue x , so the result¬ 
ing distance variation is dr = usinO = 0.15^. This is 
much smaller than the distance variation induced by the 
cohesive motion. It indicates that frictional motion (in 
the x and y-directions) results in very small variations in 
the LJ potential, which is the underlying mechanism for 
the friction issue in the interlayer LJ potential for layered 
materials. 

We find that the above friction issue for the LJ poten¬ 
tial in layered materials can be eliminated by introducing 
only one more energy parameter. We propose to supple¬ 
ment the LJ potential by the following Gaussian shaped 
potential, 

V G =ge~p\ (2) 

where p = \Jx 1 + y 2 is the projection of the distance onto 
the xy-plane, and where g is the only parameter for the 
Gaussian potential. The physical essence of this Gaus¬ 
sian potential is to guarantee the AA-stacking graphene 
layers to be the highest-energy configuration. This Gaus¬ 
sian potential impacts the frictional energy, but has no 
effect on the interlayer cohesive energy in the layered ma¬ 
terials. Fig. [5] shows the Gaussian potential curve. The 
strength of the relative frictional motion is directly re¬ 
lated to the frequency of the C mode, so the parameter 
g = 94.87 meV is determined by fitting to the frequency 
of the C mode in bilayer graphene. 

The total interaction energy between graphene layers 
is a combination of the LJ potential and the Gaussian 
potential, 

V = V LJ + Vg, (3) 

where the LJ portion is calculated in Eq. m and the 
Gaussian portion is calculated in Eq. ©• There are in 
total three potential parameters, with e and a in the LJ 
potential and g in the Gaussian potential. 


TABLE II: LJG parameters for bilayer M 0 S 2 . Numbers in 
the parentheses are the experimental lattice constant^ and 
frequency.— The first row is for the LJ potential (g = 0), while 
the second row is the LJ-G potential. Energy parameter is in 
meV. Length parameter is in A. Frequency is in cm -1 . 


e 

a 

g 

c (12.3) 

cob (40.2) 

wc (23.1) 

23.6 

3.18 

0 

12.36 

40.2 

14.5 

23.6 

3.18 

175.68 

12.36 

40.2 

23.1 


TABLE III: LJG parameters for bilayer BP. Numbers in 
the parentheses are experimental lattice constant^ and 
frequency.— The first row is for the LJ potential (g — 0), 
while the second row is the LJ-G potential. Energy param¬ 
eter is in meV. Length parameter is in A. Frequency is in 
cm -1 . 


e 

a 

g (meV) 

c (10.478) 

ujb (61.6) 

cocx (13.7) 

WCy (36.5) 

15.94 

3.438 

0 

10.5254 

59.3 

16.15 

18.11 

15.94 

3.438 

123.0 

10.5254 

59.3 

29.0 

36.2 


Following the same procedure as bilayer graphene, we 
can obtain LJ-G potential for M0S2 bilayers and BP bi¬ 
layers. Tab. HH lists LJ-G potential parameters for the 
M0S2 layers, while Tab. IIIII shows the LJ-G potential 
parameters for the BP layers. The fitting procedure to 
obtain the three parameters cr, e and g is the same as for 
graphene, i.e. they were obtained by fitting to the out-of¬ 
plane lattice constant, interlayer breathing mode, and in¬ 
terlayer shear mode, respectively. It should be noted that 
BP is highly anisotropic in the two in-plane directions 
resulting from its puckered configuration, which leads to 
different frequencies for the two C modes in bilayer BP^ 
The LJ-G potential can only provide an accurate descrip¬ 
tion for one C mode, since there is only one potential 
parameter in the Gaussian potential. Potentials with at 
least two parameters (eg. two independent Gaussian po¬ 
tentials) are needed to describe accurately both C modes 
in BP layers. 

One advantage of the LJ-G potential proposed here is 
that the potential parameters for heterostructures con¬ 
structed using different layered materials can be ex¬ 
tracted using the standard geometric combination rules 

e = V e i 6 2 


A — A 1 A 2 . 

Hence, the LJ-G potential can be easily applied to study 
the interlayer interactions in graphene/MoS 2 /BP het- 
erostructures. 


III. COHESIVE AND FRICTIONAL ENERGY 

In the previous section, we have proposed the LJ-G 
potential to describe the interlayer interaction of layered 
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FIG. 6: (Color online) Energy variations for bilayer graphene 
from the LJ and Gaussian potentials. Top panel: Gaussian 
potential has no contribution to the cohesive energy. Middle 
and bottom panels: Gaussian potential contributes 96.1% of 
the frictional energy along x and y directions. The spring in 
the inset indicates the interlayer interaction. 


materials. The rest of this paper is devoted to the appli¬ 
cation of the LJ-G potential for computing the cohesive 
energy and frictional energy in graphene, M 0 S 2 , BP, and 
their heterostructures. 

Fig. [6] top panel shows the cohesive energy in bilayer 
graphene. The zero energy point is set at the equilib¬ 
rium interlayer spacing for bilayer graphene. The co- 




FIG. 7: (Color online) Energy variations for bilayer M 0 S 2 
from LJ potential and Gaussian potential. Top panel: Gaus¬ 
sian potential has no contribution to the cohesive energy. 
Middle and bottom panels: Gaussian potential contributes 
60.3% of the frictional energy along x and y directions. The 
spring in the inset indicates the interlayer interaction. 


hesive energy is contributed solely by the LJ potential, 
while the Gaussian potential has no contribution. The 
cohesive energy in the limit of large interlayer spacing is 
24.3 meV/atom, which can be regarded as the cohesion 
energy for bilayer graphene. This cohesive energy value 
is in good agreement with the value of 23.0 meV/atom 
from the first principles calculations^ The frictional en- 
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FIG. 8: (Color online) Energy variations for bilayer BP from 
LJ potential and Gaussian potential. Top panel: Gaussian 
potential has no contribution to the cohesive energy. Middle 
and bottom panels: Gaussian potential contributes 75.0% of 
the frictional energy along x and y directions. The spring in 
the inset indicates the interlayer interaction. 


FIG. 9: (Color online) Energy variations for graphene/MoS 2 , 
graphene/BP, and M 0 S 2 /BP heterostructures. The interlayer 
interaction is described by the LJ-G potential. Top panel: co¬ 
hesive energies for the heterostructures are almost the same as 
the individual compositions. Middle and bottom panels: fric¬ 
tional energies for the heterostructures are one order smaller 
than the individual constituents. The spring in the inset in¬ 
dicates the interlayer interaction. 


ergy curves along the x and y directions are shown in 
the middle and bottom panels. It can be seen that 
the LJ potential makes a very limited contribution to 
the total frictional energy, which is dominated by the 
Gaussian potential. From the middle panel, the AA- 
stacking bilayer graphene has the maximum frictional 
energy of 15.0 meV/atom, which is in the range of pre¬ 


viously reported values of 13.0 meV/atom in Ref@ and 
19.0 meV/atom in Ref [5|. 

Fig. [7] shows the cohesive energy and the frictional en¬ 
ergy in bilayer M 0 S 2 . The cohesive energy for bilayer 
M 0 S 2 is 22.6 meV/atom, which is quite close to the cohe- 
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sive energy of bilayer graphene. For the frictional energy, 
an obvious difference between bilayer M 0 S 2 and bilayer 
graphene is that the contribution from the LJ potential 
for the frictional energy in bilayer M 0 S 2 is 39.7%, which is 
considerably larger than 3.9% in bilayer graphene. The 
contribution percentage of each individual potential is 
computed based on the energy variation around the en¬ 
ergy minimum point. Furthermore, the maximum fric¬ 
tional energy for bilayer M 0 S 2 is significantly larger than 
bilayer graphene, reaching about 40 meV/atom in the 
x-direction and about 15 meV/atom in the y-direction. 

The cohesive energy and frictional energy for the bi¬ 
layer BP is shown in Fig. [U The cohesive energy 
is 32.1 meV/atom which is larger than both bilayer 
graphene and M 0 S 2 . The LJ potential contributes 
25% to the frictional energy in bilayer BP. The maxi¬ 
mum frictional energy for bilayer BP is in between that 
of bilayer graphene and bilayer M 0 S 2 , reaching about 
33 meV/atom in the x-direction and about 22 meV/atom 
in the y-direction. 

For the graphene/MoS 2 , graphene/BP, and M 0 S 2 /BP 
heterostructures, the LJ-G potential parameters are de¬ 
termined by the combination rule in Eq (|5j). Fig. [9] top 
panel shows that the cohesive energy of all three het- 
erostructures are very similar. The frictional energy for 
these heterostructures are at least one order smaller than 
each individual constituent. The weak frictional energy 
in the heterostructure is due to the intrinsic lattice mis¬ 
match of the two individual constituents^ This weak 
frictional energy can also be analyzed from a geometrical 
point of view^ The intrinsic lattice mismatch leads to 
a Moire pattern, resulting in a large unit cell for the 
heterostructure. The Moire pattern varies during the 
frictional motion of the heterostructure. The large unit 


cell contains lots of inequivalent atoms; i.e., these atoms 
have different contribution to the interlayer interaction. 
The total interlayer energy is the summation over the 
potential for all of these inequivalent atoms. Mathemat¬ 
ically, the summation can be regarded as an integration, 
which is independent of the details for the Moire pattern. 
Hence, the interlayer potential remains almost unchanged 
during the frictional motion. 


IV. CONCLUSION 

In conclusion, we have proposed a Gaussian potential 
with only one parameter to supplement the standard LJ 
potential in layered materials such as graphene, M 0 S 2 , 
BP, and their heterostructures. The Gaussian potential 
governs the frictional motion of the layered system, while 
it has no effect on the cohesive motion for layered ma¬ 
terials. The LJ-G potential energy parameters are fitted 
to the frequency of the interlayer B mode and C mode. 
As an application of the LJ-G potential, we calculated 
the interlayer cohesive energy and frictional energy in 
graphene, M 0 S 2 , BP, and their heterostructures. Due 
to the intrinsic lattice mismatch in the heterostructure, 
the frictional energy for the heterostructure is found to 
be one order smaller than the frictional energy in each 
individual constituent. 
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